function hfile = srf_BackToVMR(hfile, opts)
% SRF::BackToVMR  - back-project vertices to VMR space
%
% FORMAT:       [srf = ] srf.BackToVMR([opts])
%
% Input fields:
%
%       opts        optional settings
%        .fillmode  either of {'nearest'} or 'linear'
%        .from      from V + from * N (default: -0.75)
%        .step      step (default: 0.5)
%        .to        to V + to * N (default: 0.75)
%
% Output fields:
%
%       srf         SRF object with VertexVMRData set accordingly

% Version:  v0.7g
% Build:    9040304
% Date:     Apr-03 2009, 4:00 AM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
if nargin < 2 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'fillmode') || ...
   ~ischar(opts.fillmode) || ...
   ~any(strcmpi(opts.fillmode(:)', {'linear', 'nearest'}))
    opts.fillmode = 'nearest';
else
    opts.fillmode = 'linear';
end
if ~isfield(opts, 'from') || ...
   ~isa(opts.from, 'double') || ...
    numel(opts.from) ~= 1 || ...
    isinf(opts.from) || ...
    isnan(opts.from) || ...
    opts.from < -8 || ...
    opts.from > 8
    opts.from = -0.75;
end
if ~isfield(opts, 'step') || ...
   ~isa(opts.step, 'double') || ...
    numel(opts.step) ~= 1 || ...
    isinf(opts.step) || ...
    isnan(opts.step) || ...
    opts.step < 0.125 || ...
    opts.step > 4
    opts.step = 0.5;
end
if ~isfield(opts, 'to') || ...
   ~isa(opts.to, 'double') || ...
    numel(opts.to) ~= 1 || ...
    isinf(opts.to) || ...
    isnan(opts.to) || ...
    opts.to < opts.from || ...
    opts.to > 8
    opts.to = max(opts.from, 0.75);
end

% get content
bc = bvqxfile_getcont(hfile.L);
c = bc.VertexCoordinate;
n = bc.VertexNormal;

% create VMRData
vmrd = uint8([]);
vmrd(256, 256, 256) = uint8(0);

% iterate over range
for sc = opts.from:opts.step:opts.to
    
    % for nearest
    if strcmpi(opts.fillmode, 'nearest')
        
        % get coordinate (+ sc * n)
        tc = round(c + sc * n);
        
        % remove out of bounds
        tc(any(tc < 1 | tc > 256, 2), :) = [];
        
        % get index
        tx = sub2ind(size(vmrd), tc(:, 1), tc(:, 2), tc(:, 3));
        
        % and set to 235
        vmrd(tx) = 235;
        
    % for linear
    else
        
        % get coordinate
        tc = c + sc * n;
        
        % get lower bound and diff
        tcl = floor(tc);
        tc = tc - tcl;
        
        % remove all that cannot make it
        tcx = any(tcl < 1 | tcl > 255, 2);
        tcl(tcx, :) = [];
        tc(tcx, :) = [];
        
        % iterate for x, y, z
        for xd = 0:1
            for yd = 0:1
                for zd = 0:1
                    
                    % get index
                    tcx = sub2ind(size(vmrd), ...
                        tcl(:, 1) + xd, tcl(:, 2) + yd, tcl(:, 3) + zd);
                    
                    % set to higher value
                    vmrd(tcx) = max(vmrd(tcx), uint8(25 + round(200 .* prod(...
                        [abs(tc(:, 1) - xd), abs(tc(:, 2) - yd), abs(tc(:, 3) - zd)], 2))));
                end
            end
        end
    end
end

% put back into VOI and return!
bc.VertexVMRData = vmrd;
bvqxfile_setcont(hfile.L, bc);
